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We perform lattice simulations of two flavors QCD using the optimal domain-wall fermion, in 
which the chiral symmetry is preserved to a good precision (nires ^ 0.3 MeV) on the 16^ x 32 
lattice (L ^ 2 fm) with inverse lattice spacing a^^ ^ 1.8 GeV, and N^ — 16 in the fifth dimension, 
for eight sea quark masses corresponding to the pion masses in the range 210-500 MeV. We 
present our first results of the mass and the decay constant of the pseudoscalar meson, which are 
in good agreement with the next-to-leading order chiral perturbation theory for Mj^ < 450 MeV, 
and from which we determine the low-energy constants /, E, Ij and I4. At the physical pion 
mass Mji = 135 MeV, we obtain the pion decay constant fj^ = 133(1)(2) MeV, and the average 
up and down quark mass m^f{2 GeV) = 4.09(7) (11) MeV, where the first error is statistical, 
and the second error is systematic due to the truncation of the higher order corrections and the 
uncertainty in the determination of the lattice spacing. Furthermore, we also obtain the chiral 
condensate 2:^(2 GeV) = [250(4)(7)MeV]l 
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1. Introduction 

Lattice QCD with exact chiral symmetry is an ideal theoretical framework to study the non- 
perturbative physics from the first principles of QCD. However, it is rather nontrivial to perform 
Monte Carlo simulation such that the chiral symmetry is perserved to a very high precision and all 
topological sectors are sampled ergodically. 

Since 2009, the Taiwan Lattice QCD Collaboration (TWQCD) has been using a GPU cluster 
(currently constituting of 250 NVIDIA CPUs) which attains 40 Teraflops (sustained) to simulate 
unquenched lattice QCD with the optimal domain- wall quarks [jl|, We have realized our goal of 
preserving the chiral symmetry to a good precision (with m^s ~ 0.3 MeV) and also sampling all 
topological sectors ergodically. 

In this paper, we present our first results of the mass and the decay constant of the pseudoscalar 
meson in two flavors QCD, and compare our results with the next-to-leading order (NLO) chiral 
perturbation theory (ChPT). We find that our data is in good agreement with NLO ChPT for M,^ 
less than 450 MeV, and from which we determine the low-energy constants /, £, 1^ and T4, and 
and the average up and down quark mass m^f{2 GeV). Our result of the topological susceptibility 
is presented in Ref. ||3|], and our strategy of using GPU to speed up our Hybrid Monte Carlo 
simulations is presented in Ref. [Q]. 

2. Hybrid Monte Carlo Simulation with Optimal Domain-Wall Quarks 

The optimal domain-wall fermion is the theoretical framework which preserves the (mathe- 
matically) maximal chiral symmetry for any finite Ng (the length of the fifth dimension). Thus the 
artifacts due to the chiral symmetry breaking with finite A^v can be reduced to the minimum. 

The action of the optimal domain-wall fermion is defined as [|l|] 

5odwf= £'/^r.v[(w,Dw + l)«'5,y + (C0,D„,-l);^/L,y]v/yy (2.1) 

where the weights {(0^} along the fifth dimension are fixed according to the formula derived in Ref. 
^ such that the maximal chiral symmetry is attained. Here D^. denotes the standard Wilson-Dirac 
operator plus a negative parameter —mo (0 < mo < 2), 

(£>w)x.' = -^£[(l-7M)^MW4+Ay + (l + rM)f^M(^')4-A/] +(4-'«o), (2.2) 



2V 



and 



L = P+L++P L , P± = (l±75)/2, (2.3) 



-{mq/2mo)ON^^,f, s= 1 

where m^ denotes the bare quark mass. Separating the even and the odd sites on the 4D space-time 
lattice, ^ can be written as 
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where 

Ms = 



(4-mo) + V« \l-L){l+L) ^ , {(o),s' = (oA,', (2.6) 

Si=M5Vco'\ S2 = {l+L)-^Vco'\ (2.7) 



and the Schur decomposition has been used in the last equality of (2.5), with the Schur compliment 



C=l-M5D°EM5Df°. (2.8) 

Since det^ = det5j"' • detC • det52 and and ^2 do not depend on the gauge field, we can 
just use C for the HMC simulation. After including the Pauli-Villars fields (with niq = 2mo), the 
pseudo-fermion action for 2-flavor QCD (m„ = m^) can be written as 

Spf = fclyiCC'^y^Cpv^, Cpv = C(2mo). (2.9) 

In the HMC simulation, we first generate random noise vector ^ with Gaussian distribution, 
then we obtain (j) = CpyCt, using the conjugate gradient (CG). With fixed 0, the system is evolved 
with a fictituous Hamiltonian dynamics, the so-called molecular dynamics (MD). In the MD, we 
use the Omelyan integrator [^], and the Sexton-Weingarten multiple-time scale method ||^. The 
most time-consuming part in the MD is to compute the vector Tj = (CC^)^'Cpi/0 with CG, which is 
required for the evaluation of the fermion force in the equation of motion of the conjugate momen- 
tum of the gauge field. Thus we program the GPU to compute tj , using CG with mixed precision 
[0]. Also, we have ported the computation of the gauge force and the update of the gauge field to 
the GPU. 

Furthermore, we introduce an extra heavy fermion field with mass itie {mq <^ nin < 2mo), 
similar to the case of the Wilson fermion [0]. For two flavors QCD, the pseudofermion action 
becomes 

Spf = <^^cl--^CH<p + (^lcly-^Cpv(pH, CH = C{mH), (2.10) 



which gives exactly the same fermion determinant of (2.9). However, the presence of the heavy 
fermion field plays the crucial role in reducing the light fermion force and its fluctuation, thus 
diminishes the change of the Hamiltonian in the MD trajactory, and enhances the acceptance rate. 

For a system with CPU and GPU, we can have both of them compute concurrently, e.g., while 
the GPU is working on the CG of the light quark field, the CPU can compute the fermion force 
of the heavy fermion field. This asynchronous concurrent excecution mode enhances the overall 
performance by ~ 5%. 

A detailed description of our HMC simulations will be presented in a forthcoming paper [0]. 



3. Lattice setup 

We simulate two flavors {Nf = 2) QCD on the 16-' x 32 lattice at the lattice spacing a ~ 0.11 fm, 
for eight sea quark masses in the range m^a = 0.01,0.02, • • • ,0.08. For the gluon part, we use the 
plaquette action at j8 = 5.90. For the quark part, we use the optimal domain-wall fermion with 
Ns = 16. After discarding 300 trajectories for thermalization, we accumulated about 3000 — 3200 
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Figure 1: The topological charge versus the trajectory in the HMC simulation of two flavors QCD with 
ODWF. The lattice is 16^ x 32 with the spatial box ~ (2 fm)^, and the quark mass corresponding to M;^ ^ 300 
MeV. The topological charge is obtained by projecting the zero modes of the overlap Dirac operator. 

trajectories in total for each sea quark mass. From the saturation of the error (by binning) of the 
plaquette, as well as the evolution of the topological charge (see Fig. |l]), we estimate the autocorre- 
lation time to be ~ 10 trajectories. Thus we sample one configuration every 10 trajectories. Then 
we have 270 — 290 configurations for each sea quark mass. 

We determine the lattice spacing by heavy quark potential with Sommer parameter = 0.49 
fm. The inverse lattice spacing versus the quark mass is plotted in Fig. ^ Using the hnear fit, we 
obtain the inverse lattice spacing in the chiral limit, = 1.8153(28) GeV. 
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Figure 2: The inverse lattice spacing a [GeV] versus niqa for two flavors QCD with ODWF. 

For each configuration, we calculate the exact zero modes plus 80 conjugate pairs of the 
lowest-lying eignmodes of the overlap Dirac operator. We outline our procedures as follows. First, 
we project 240 low-lying eigenmodes of using v-TRLan alogorithm [P], where each eigenmode 
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has a residual less than 10^^^. Then we approximate the sign function of the overlap operator by 
the Zolotarev optimal rational approximation with 64 poles, where the coefficents are fixed with 
^max — (6-4)^, and equal to the maximum of the 240 projected eigenvalues of H^^,. Then the 
sign function error is less than 10^^^. Using the 240 low-modes of H}^, and the Zolotarev approxi- 
mation with 64 poles, we project the zero modes plus 80 conjugate pairs of the lowest-lying eign- 
modes of the overlap operator with the v-TRLan algorithm, where each eigenmode has a residual 
less than 10"^^. 

We measure the chiral symmetry breaking (due to finite A'^^) by computing the residual mass 

{[/} \ ^[{dI + m^) {Dc + m^jJo / 



where (D^ + m^)^^ is the valence quark propagator with ruq equal to the mass of the sea quark, tr 
denotes the trace running over the color and Dirac indices, and the subscript {U} denotes averaging 
over an ensemble of gauge configurations. It turns out that, after averaging over an ensemble of a 
few hundreds of independent gauge configurations, nires is insensitive to the location of the origin 



= (0,0,0,0). Thus (3.1) gives a reliable measure of chiral symmetry breaking due to finite A^^. 



The derivation of ( |3.1| ) will be given in a forthcoming paper [|10|]. 

In Fig. ^, we plot the residual mass versus the quark mass. Using the power-law fit, we obtain 
the residual mass in the chiral limit, niresa = 0.00018(2), which amounts to rrires = 0.32(4) MeV. 
Note that the value of rures is less than 1/10 of the statistical and systematic errors of the inverse 
lattice spacing, thus confirming that the chiral symmetry has been preserved to a good precision in 
our simulation. 

10-3 1 1 




Figure 3: The residual mass versus the quark mass for two flavors QCD with ODWF. 



4. The Mass and the Decay Constant of the Pseudoscalar Meson 

In this section, we present our first results of the pseudoscalar mass and decay constant, for 
2 flavors QCD with optimal domain- wall quarks and the plaquette gluon action at j8 = 5.90, on 
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the 16^ X 32 X 16 lattice. In Fig. ^ we plot M^/mq and /;;: versus nig respectively. Here we have 
made the correction for the finite volume effect using the estimate within ChPT calculated up to 
^{M'^/{4nfj[)^) [12], since our simulation is done on a finite volume lattice with Mj^L ~ 2.0 for 



the lightest sea quark, and its finite volume effect cannot be neglected. 

Taking into account of the correlation between M^/mg and fji for the same sea quark mass, we 
fit our data to the formulas of the next-to-leading order (NLO) chiral perturbation theory (ChPT) 



[11] 
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where A, are related to the low energy constants 



h = In 



L = In 



0.140 GeV. 



(4.1) 
(4.2) 

(4.3) 



For the six lightest quark masses (corresponding to pion masses in the range 210 — 445 MeV), 
our fit gives 

1 = 0.2105(30) GeV, / = 0.127(2) GeV, Z3= 4.37(18), [4 = 5.31(11), (4.4) 

with ;i;^/dof = 0.4. At the physical pion mass M„ ~ 0.135 GeV, the value of pion decay con- 
stant is fji = 0.133(1) GeV, and the bare quark mass is 0.0069(2) GeV. In order to convert the 
bare quark mass to that in the MS scheme, we calculate the renormalization factor Z^^ (2 GeV) 



using the non-perturbative renormalization technique through the RI/MOM scheme []13[], and ob- 
tain Z^^{2 GeV) = 0.5934(10) [p^. Then the value of the average up and down quark mass is 
transcribed to 



m^f{2 GeV) =4.09(7) (11) MeV, 
Similarly, the value of L in ( p^ is transcribed to 

L^{2 GeV) = [250(4) (7) MeV]^ 



(4.5) 



(4.6) 



The systematic error is estimated from the turncation of higher order effects and the uncertainty in 
the determination of lattice spacing with ro = 0.49 fm. Since our calculation is done at a single 
lattice spacing, the discretization error cannot be quantified reliably, but we do not expect much 
larger error because our lattice action is free from 0{a) discretization effects. 



5. Concluding remarks 

Using a GPU cluster (currently attaining 40 sustained Teraflops with 250 NVIDIA GPUs), 
we have succeeded to simulate unquenched lattice QCD with optimal domain-wall quarks, which 
preserves the chiral symmetry to a good precision and samples all topological sectors ergodically. 
Our results of the mass and the decay constant of the pseudoscalar meson (in this paper) and the 
topological susceptibility (in Ref. [^) suggest that the nonperturbative chiral dynamics of the sea 
quarks are well under control in our simulations. This provides a new strategy to tackle QCD 
nonperturbatively from the first principles. 
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Figure 4: Physical results of 2-flavor QCD with optimal domain-wall quarks: (a) m^/m^, and (b) fj^. The 
solid lines are the simultaneous fits to the NLO ChPT, for the six lightest quark masses. 
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